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Abstract. Let / be a univariate polynomial with real coefficients, / £ R[X], 
Subdivision algorithms based on algebraic techniques (e.g., Sturm or Descartes 
methods) are widely used for isolating the real roots of / in a given interval. 
In this paper, we consider a simple subdivision algorithm whose primitives are 
purely numerical (e.g., function evaluation). The complexity of this algorithm 
is adaptive because the algorithm makes decisions based on local data. The 
complexity analysis of adaptive algorithms (and this algorithm in particular) 
is a new challenge for computer science. In this paper, we compute the size of 
the subdivision tree for the SqFreeEVAL algorithm. 

The SqFreeEVAL algorithm is an evaluation-based numerical algorithm which 
is well-known in several communities. The algorithm itself is simple, but prior 
attempts to compute its complexity have proven to be quite technical and 
have yielded sub-optimal results. Our main result is a simple 0(d(L + lnd)) 
bound on the size of the subdivision tree for the SqFreeEVAL algorithm on the 
benchmark problem of isolating all real roots of an integer polynomial / of 
degree d and whose coefficients can be written with at most L bits. 

Our proof uses two amortization-based techniques: First, we use the alge- 
braic amortization technique of the standard Mahler-Davenport root bounds 
to interpret the integral in terms of d and L. Second, we use a continuous 
amortization technique based on an integral to bound the size of the sub- 
division tree. This paper is the first to use the novel analysis technique of 
continuous amortization to derive state of the art complexity bounds. 
Key words: Continuous Amortization, Adaptive Analysis, Subdivision Algo- 
rithm, Integral Analysis, Amortization, Root Isolation. 



1. Introduction 

In this paper, we show that the size of the subdivision tree for the simple, 
evaluation-based, numerical algorithm SqFreeEVAL has size 0(d(L + hid)) for the 
benchmark problem of isolating all of the real roots of an integer polynomial of 
degree d whose coefficients can be represented by at most L bits. Under the mild 
assumption that L > hid, this complexity simplifies to the optimal size of O(dL). 
The optimality and simplicity of the SqFreeEVAL algorithm imply that it may be 
a useful algorithm in practical settings. The bound on the size of the subdivision 
tree is achieved via a straight-forward and elementary argument. The two main 
techniques which are used in the computation are algebraic amortization, in the 
form of Mahler-Davenport bounds, and continuous amortization, in the form of an 
integral technique as presented in 

1.1. EVAL-type algorithms. The SqFreeEVAL algorithm which we study in this 
paper is a specific example of what we call an EVAL-type algorithm. These algo- 
rithms are so named because they are based on function evaluation: EVAL-type 
algorithms take, as input, functions which allow some subset of the following two 
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predicates: First, these functions and their derivatives can be evaluated at a count- 
able dense subset of their domain. In this paper, the domain will be the real 
numbers and the countable dense subset will be the dyadic integers. Second, these 
functions and their derivatives can be approximated on intervals in such a way that 
the approximation converges as the input intervals converge to a point. In this 
paper, the approximation is derived from interval arithmetic on a Taylor sequence. 
The simplest and most well-known ex ample of an EVAL-type algo rithm is Lorensen 
and Cline's marching cube algorithm Lorensen and Clind . Il987| . 

EVAL-type algorithms are typically studied because of their simplicity and gen- 
erality. These algorithms are fairly general because their inputs can be extended 
to more general analytic functions. In particular, many analytic functions have 
interval arithmetic available to them, and, therefore, it is possible to approximate 
these functions on intervals. In addition, with the limited predicates available to 
EVAL-type algorithms, most of the techniques which are used in these algorithms 
are analytically based (as opposed to algebraically based). These algorithms are 
simple because, in many cases, EVAL-type algorithms are based on simple recur- 
sive bisection algorithms. Such algorithms iteratively subdivide an initial domain 
until each set in the resulting partition of the initial domain satisfies a (usually 
simple) terminal condition . Bisection algorithms are common in computer graphics 
Boier-Martin et all 20051 as well as in computational science and engine ering appli- 
cations International Conference on Domain Decomposition Methods! . Bisection 
algorithms are of particular interest because they are adaptive; they perform more 
bisections near difficult features and fewer bisections elsewhere. However, this adap- 
tivity makes the complexity analysis of such algorithms more difficult because the 
subdivision tree may have a few deep paths while the remainder of the tree remains 
modest in size. 
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ate case in [Galehousd . [2009] . All of these algorithms are devoted to approximating 
algebraic (an d in some cases analytic varieties) in the real or complex settings. The 
algorithms in [Burr et ah , 2011 , 20091 ] are designe d to find all real roots of a polyno- 
mial or analytic function while the algorithms in [Henrici Il970l lYakoubsohnL 120051 
ISagraloff and Yarl |2009| are designed to fin d the complex roots of a polynomial or 
analytic function (note that [Henrici , 197Cl| is only designed to find a single root of 
a polynomial). Each of these algorithms is very closely related to the SqFreeEVAL 
algorithm considered in this paper; the main differences are in the setting, in the 
type of subdivisions performed, and in various preprocessing steps. We give a more 
detailed account of t hese algorithms in the next section. The t wo-dimensional 
EVAL-type algorithm Plantinga and Vegter . 2004 . Plantinga , 2006j was presented 
for approximating smoo th and bounded v arieties. It was extended to singular and 
unbounded varieties in IBurr et all 120101 : in ad dition, the tests performed by the 
algorithm were improved in |Lin and Yap . l2009l j . 



1.2. The SqFreeEVAL algorithm. There are many bisection algorithms for finding 
roots, see Section ITTSI for references, but among such algorithm s, the Sq FreeEVAL 



algorithm is one of the simplest and most widely applicable, see [Burr et all 12011 1 



There are two distinct paths in the literature which arrive at algorithms similar 
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to the SqFreeEVAL algorithm: one path proceeds through the consideration of 
magnitudes of derivatives and the other path proceeds via interval arithmetic. 
We b egin by discuss ing the history from the magnitudes of derivatives perspec- 



tive. In Henrici 



1970j |. ts an algorithm for finding a single com- 

plex root of a poly nomial. The test T3 from the paper is essentially used here. In 
[Yakoubsohnll2005j . the test is developed into a bisection algorithm and to find all 
complex r oots of entire f unctions, not just polynomials. In the paper, however, the 
test from Henrici . 19TCl| is used only as a one-sided test; therefore, the algorithm 
can only exclude regions from containing roots and does not confirm that roots 
exist in the final regions. There, the algorithm was termed a bisection- exclusion 
algorithm to reflect this drawbac k. Finally, in |Sagraloff and Yap . l2009l |. the al- 
gorithm from [Yakoubsohn , 2005j | was adapted to polynomials in order to confirm 
that roots exist in the final regions; there, the authors studied both an algorithm 
for finding complex roots as well as one for finding real roots. The SqFreeEVAL 
algorithm is a natural restriction of these complex root-finding algorithms to the 
real line. 

On the other hand, from the interval arithm etic community, a bisection algo- 
rithm using interval methods was suggested in [MooreL Il96fil iMitchehl . flQQnj | . In 
these papers, any interv al function can be used; if t he standard centered form for 
polynomials is used, see Ratschek and Roknel . 1984 1 . then the exclusion conditions 
are identical (when / and /' are square free) to those for the SqFreeEVAL algorithm. 

In this paper, we study the SqFreeEVAL algorithm on the standard benchmark 
problem of finding all of the real roots of a polynomial. We show that, in this 
case, the subdivision tree has the favorable size of 0(d(L + Ind)) which simplifies 
to the optimal size of 0(dL) under the mild assumption that L > hid. Since this 
algorithm uses only local data to find roots, it is an adaptive algorithm and may be 



more efficient than the standard exact algorithms in certain cases, see [Burr et al 
120111. In additio n, the SqFreeEVAL algorithm can handle analytic varieties, see 



[Burr et all 12011 1 , which extends its reach beyond that of more standard exact 



algorithms which require sophisticated algebraic primitives and are specialized to 
polynomials. These advantages of the SqFreeEVAL algorithm imply that it may be 
more practical than other standard root isolation algorithms in practice. 



1.3. Previous complexity results. The computational complexity of EVAL-type 
algorithms has proven to be quite a challenging problem because the algorithms 
are adaptive and the analytic primitives do not carry much information about the 
global structure (unlike algebraic information). Here, we survey the complexity 
analyses of the precursors to the SqFreeEVAL algorithm. In most situations, the 
complexity is computed in terms of the size of the subdivision tree of the specific 
EVAL-type algorithm (this is almost equivalent to counting the number of tests 
performed by the algorithm). There have been two main techniques to find the 
size of the subdivision tree: by finding the width of the subdivision tree at various 
subdi vision levels or by finding the local depth of the subdivision tree. 

In [Henrici Il970| . the author is searching for only a single root, and, therefore, 
retains a single disk containing a root at each stage of the algorithm. Many tests 
are performed in the algorithm, however, because at each stage of the algorithm, 
tests are performed on a covering of the previously retained disk. The final stopping 
criterion for this algorithm is based on a precision e > which is chosen a priori 
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by the user. When the worst-case root separation bound for a polynomial is used, 
the c omplexity of the sub division tree becomes 0(d 3 (L + lnc?)). 

In Yakoubsohn , 2005], the author is searching for all of the complex roots of an 



analytic function. In the computation, a bound on the width of the tree is computed 
to bound the number of subdivisions performed. Since this algorithm only excludes 
regions and lacks an inclusion test, it is possible that the final output regions do 
not contain roots or do not separate roots. The final stopping criterion for this 
algorithm is based on a precision e > which is chosen a priori by the user. When 
the worst-case root separation bound for a polynomial is used, the complexity of 
the subdivision tree becomes either 0(d 4 (L + lnc?)) or 0(d 3 (L + In 3 d)) after [In of] 
steps of the Graeffe iter ation. 

In Burr et al. . 2009j |. we search for all of the real roots of a polynomial. Here, 



the computation is based on the depth of the tree over each point of the initial 
interval. In the paper, we introduced the idea of continuous amortization via an 
integral and showed how to use it to bound the size of the subdivision tree. In 
particular, we proved a complexity bound of 0(d 3 (L + hid)) for the subdivision 
tree. 

In Sagraloff and Yapl . l2009j | , the authors present algorithms to find all of the 



real or all of the complex roots of a polynomial. In the computation, a bound 
on the width of the subdivision tree is used to compute the number of subdivi- 
sions performed. The authors show that the complexity of the subdivision tree is 
0(d(L + lnd)(lnL + hid)) in the real case and 0((d\nd) 2 (L + \nd)) in the com- 
plex case. In addition, the authors show that the bit complexity of both of these 
algorithms is 0(d 4 L) where the O means that logarithmic factors in d and L have 
been suppressed. 

Each of the analyses in Yakoubsohn . 2005 . Burr et al. . 20091 Sagraloff and YaS 



|2009| are quite technical, complicated, and require several constants to be defined 
whose use becomes justified only after the completion of the complexity analysis. 
In contrast, the computation in this paper is quite simple and provides the better 
bound of 0(d(L + \nd)). It should be noted that although this is the best bound 
known, it does not directly replace the bounds presented in these papers because 
some are in the different setting of the complex plane and others use different 
preprocessing steps. In the case where the polynomial and its derivative are both 
square free and we are searching for the real roots, all of these algorithms are 
identical and our bound on the subdivision tree is the best. 

1.4. Algebraic and continuous amortization. In this paper, we use amortiza- 
tion in two forms: algebraic a nd continuous. Algebraic amortization originated with 
Davenport Davenport . 1985| where the individual root separation bou nds are re- 



placed by a product of root se parations. This bound was then studied in [Du et al 



20071 lEigenwillig et all 120061 ] where it was generalized to other root separation 
products including complex roots. This technique has proven useful to compute 
the complexity of the subdivision tree for many other root isolation techniqu es, 



see Section [T31 We introduced continuous amortization in [Burr et aL . 2009j to 



bound the size of the subdivision tree of an EVAL-type algorithm. In this paper, we 
show that continuous amortization can be used to significantly simplify complexity 
calculations. 

In continuous amortization, we use a complexity charge (f> whose domain is the 
input region, and, for each a; in the input region, 4>{x) is a lower bound on the 
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size of any leaf interval containing x. Then l/<f>(x) is related to the dep th of the 
subdivision tree for an interval which contains x. In 



Burr ct al., 2009], we used 



continuous amortization to compute the size of a subdivision tree for an EVAL- 
type algorithm. In this paper, we greatly simplify the computation and provide a 
complexity bound for the SqFreeEVAL algorithm. 

We call the function <fi mentioned in the previous paragraph a s topping function 
for the algorithm. Similar functions also appeared in lHenricil.ll970l where they were 
called inner and outer convergence functions. In Yakoubsohnf 2005] such functions 
were also termed exclusion functions. In both cases, these stopping functions were 
used to compute the complexity of the algorithm, but they were not used in a 
continuous amortization computation. 



1.5. Other root isolation algorithm s. There is an ex tensive amount of literature 
on the complexity of root isolation, see [panl . ll997l l 1996] for surveys of the previous 
literature, which we will not attempt to cover here. Most algorithms are compared 
by their performance on the benchmark problem of finding all real roots of a polyno- 
mial of degree d and whose coefficients can be represented by at most L bits. For this 
problem, the b it-complexity of Q (d 3 (L + \nd)) for complex roots was first achieved 



by Schonhage Schonhage . 1982]. In many algorithms, the size of the subdivision 



tree is smaller than this bound because, for each no de in the subdivis ion tree, ad- 
ditional calculations must be performed. Davenport (Davenportl Il985| p roved that 
the the subdivision tree for t he Sturm meth od is 0(d(L + hid)), see (Reischert . 
1997 L Lickteig and Rovl 2001 , Du et al. , 2007} . More recently, it has been shown in 



Eieenwillig et al. L 20061 that the Descartes method also achieves this bound, see 



Collins and Akritasl . Il976l lEigenwillig et all . 120061 iKrandick and Mehlhornl . 120061. 



Collins et all . 12002], These methods are optimal under the weak assumption that 
L > hid. In addition, related exact techniques using continued fractions were 
shown to have a tree size of O(dL) w hen an ideal ro ot bound is used and 0(d 2 L) 
when a more practical bound is used Sharmal l2008j . In the algebraic computing 



community, the Descarte s method appears to be one of the more practica l algo- 
rithms, see Collins et all |2002[ [johnsonl Il998i iRouillier and Zimmermann , 2004 , 
Mourrain et al. . 2005, Rou illier and Zimmermannl . 2004]. In this paper, we show 
that the subdivision tree for the SqFreeEVAL algorithm also achieves this bound; 
therefore, the SqFreeEVAL algorithm should also be considered on equal footing 
with the other more well-known root finding algorithms via the Sturm or Descartes 
methods. The SqFreeEVAL algorithm may, in addition, be considered practical 
because its computations are numerical and hence easy to implement and its sub- 
division tree has a favorable size. 



1.6. Organization of this paper. In Scction[2] we introduce the SqFreeEVAL al- 
gorithm and discuss the main condition we will use for an interval to be SqFreeEVAL 
terminal. In Section [3j we review the use of stopping functions to bound the size 
of the subdivision from [Burr et all . 2009] and create a stopping function for the 
SqFreeEVAL algorithm. In Section 0] we compute the size of the SqFreeEVAL al- 
gorithm's subdivision tree using continuous amortization via the stopping function 
technique and achieve the main result of this paper, the 0(d(L + hid)) bound on 
the size of the subdivision tree for the SqFreeEVAL algorithm. Finally, we conclude 
in Section [5] 
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2. The SqFreeEVAL algorithm 

Given an interval / = [a, b] with integer endpoints and a polynomial / with 
integer coefficients, i.e., / € Z[X], the SqFreeEVAL algorithm returns a collection of 
intervals which cover and isolate the real roots of / in (a, &), i.e., every root appears 
in an output interval and each output interval contains exactly one root (ignoring 
multiplicities). In the SqFreeEVAL algorithm, if the interval [c, d] is output, then 
(c, d) contains exactly one root of / and if [c, c] is output, then c is a root of /. 
The SqFreeEVAL algorithm maintains a (finite) partition P of the interval /, i.e., a 
finite collection of intervals whose interiors are disjoint and whose union is /. The 
SqFreeEVAL algorithm iteratively bisects the elements of P until the intervals of the 
partition P are each small enough to pass the SqFreeEVAL termination conditions 
(see Section l2~Tj) . Of interest to us is the size #P of the partition, i.e., the number 
of intervals in P. 

We begin with some terminology: For an interval J = [c, d] the width of J is 
w(J) = d — c and the midpoint of J is m(J) = (c + d)/2. Also, to bisect an 
element of the partition P means to replace the interval J = [c, d] £ P by the two 
subintervals [c, m(J)] and [m(J),d\. Note that this implies that #P is one more 
than the number of bisections done by the SqFreeEVAL algorithm, i.e., the size of 
the subdivision tree. All of the calculations done by the SqFreeEVAL algorithm will 
be performed on the dyadic integers Z[l/2] so that all of the standard operations 
are exact. This prevents well-known implementation errors from arising in practice. 



2.1. Statement of the SqFreeEVAL algorithm. In the SqFreeEVAL algorithm, 
we first replace / by its square free component, which we briefly call g. Then, 
we replace /' by its square free and relatively prime to / component, i.e., we first 
take the square free component of /' and then take the portion of this polynomial 
which is relatively prime to /. We briefly call this h. Note that g\f and h\f, and, 
moreover, the roots of g are separated by roots of h by Rolle's theorem. In the case 
where / is square free, the zeros of h partition / into monotonic regions; in the case 
where / is not square free, the zeros of h no longer have this property, but they still 
partition the roots of / (and hence the roots of g). Throughout the remainder of this 
paper, except for a brief note in Section [4.21 we use these square free substitutions 
for / and /' without mention. The bounds on the subdivision tree, however, will 
be in terms of the data for original the / and not for any replacements. 

The SqFreeEVAL algorithm creates a partition of / and determines which inter- 
vals in the partition contain roots. Initially, the partition of I is P = {I}, the 
trivial partition. 



SqFreeEVAL: AN (ALMOST) OPTIMAL REAL-ROOT ISOLATION ALGORITHM 



7 



Algorithm 2.1: The SqFreeEVAL algorithm 
Repeatedly subdivide each J € P until one of the following conditions holds 

(Co) | /(mW) |>^«m^, „ 

i=i 
d-l 

(Cx) \f'(m(J))\>J2 ,, , 

//, when subdividing, /(m(J)) = 0, then output [m(J), m(J)]. 
for eac/i interval J = [c, d] € P where Ci ZioWs and /(c) • /(d) < 0, output J 



/' \ 2 

|/^ +1 )(m(J))| /«;(J) 



The termination pro of for the SqFreeEVAL algorithm is very sim ilar to the corre 



sponding statement in [Burr et al 



qtreebvAL aigoritnm is very similar to tne corre- 
■ l2009llSagraloff and Yapl . 12009]. The correctness 



proof is slightly different from the corresponding proofs for other EVAL-typc algo- 
rithms. The correctness follows from the Taylor polynomial centered at m(J): if 
one of the conditions holds, then it follows that / (for condition Co) or /' (for 
condition C±) is never zero in J since the inequalities are equivalent to a reverse 
triangle inequality on the Taylor polynomial. The first condition implies that / has 
no zeros in J. The second condition implies that / has at most one zero in J since 
roots of /' separate zeros of / (even though / might not be monotonic due to the 
replacements above). 

2.2. SqFreeEVAL terminal intervals. In this section, we provide a sufficient con- 
dition for the SqFreeEVAL algorithm to terminate without subdividing on a given 
interval, i.e., for the interval to be SqFreeEVAL terminal. 

Definition 2.1. For any polynomial g of degree d, define a g = {a%, • • • , ad} to be 
the multiset of the roots of g. In addition, define the function E ff to be the sum of 
the reciprocals of the distances from its argument to the roots of g: 

a£a„ 



a 



Note that this function can be represented in a simple form using the harmonic 
mean HM. Then, one has 



1 



HMO 



where 



is the set of distances from x to the roots of g. 



Yif and £y will be our main objects of study. We begin with the following 
lemma which connects £/(x) and £/'(#) with conditions Co and Ci, respectively: 

Lemma 2.1. The following inequality holds for i > 0: 

/M(z) 



< [£/(*)] 



The proof is a straight-forward comput ation. See the proof of Burr et al. , 20091 
Lemma 6.2] or Sagraloff and Yap . 20091 Section 5.2] for details. 



We use this lemma to show that a simple upper bound on the width of an interval 
will ensure that the conditions in the SqFreeEVAL algorithm hold. For example, in 



8 



MICHAEL A. BURR AND FELIX KRAHMER 



condition Co, divide both sides of the inequality by \f(m(J))\ and apply Lemma 
12.11 to derive the following inequality: 

^ |/W(m(J))| f wjJl X ^ 1 { X f (m(J))w(J) 



i\\f(m(J))\ V 2 J 

If w( J) < g— j^jyj i then the sum on the RHS is bounded above by a geometric series 
with r = 1/2, and, therefore, the sum is bounded by 1. This implies that condition 
Co holds. Therefore, the condition w(J) < Y: f (m(j)) 1S sufficient to ensure that 
J is SqFreeEVAL terminal. Similarly, if w(J) < g ; (m(j)) ' then J is SqFreeEVAL 
terminal by condition Ci. 

3. Stopping functions 

In this section, we show how stopping functions can be used to compute the 
size of the subdivision tree of the SqFreeE VAL algorithm . The construction in 
Section |3"7T1 was originally presented in Burr et al. , 2009l |. but we include it here 



for completeness and because the construction in Section 13.21 requires a detailed 
understanding of the method. 

3.1. Basic properties. The use of stopping functions promises to be an important 
tool for bounding the complexity of subdivision algorithms. Most of the numerical 
algorithms appearing in the introduction may benefit from this type of analysis; 
more algorithms of this type are mentioned in the Conclusion, Section [5] We begin 
by formulating an abstract algorithm called the Bisection algorithm, which is 
intended to be the prototype of these types of algorithms in one dimension. The 
notion of stopping functions and the Bisection algorithm both easily generalize 
to higher dimensions. 

Fix a predicate B (i.e., a Boolean function) on intervals with the following prop- 
erty: if K C J and B(J) is true, then B(K) is also true. The Bisection algorithm 
is the following algorithm: given an interval /, the algorithm maintains a parti- 
tion P of /. Initially, let the partition be the trivial partition P = {1} and let 
-PBisection(-f) be the final partition. 



Algorithm 3.1: The Bisection algorithm 
Repeatedly subdivide each J £ P until the following condition holds: 
B(J) is true 



A stopping function for the Bisection algorithm with predicate B is a real- 
valued function F with the following property: if, for a given interval J, there 
exists a point p e J such tha t w(J) < F(p), then B{J) is true. The following 
theorem, which also appears as Burr et al. , 20091 Theorem 3.5], bounds the number 



of subdivisions performed by the Bisection algorithm. 

Theorem 3.1. Burr et all 120091 Theorem 3.5] Let F be a stopping function for 



the Bisection algorithm, then 



( f 2dx } 

#-PBisection(j) < max j 1 , J -^-y j> 



If the Bisection algorithm does not terminate, then the integral is infinite. 



SqFreeEVAL: AN (ALMOST) OPTIMAL REAL-ROOT ISOLATION ALGORITHM 



9 



Proof. If #PBisection = 1, then the bound is immediate. If #PBisection > 1, then 
an examination of the Bisection algorithm shows that for J G ^Bisection there is 
a lower bound on w(J) since the Bisection did not terminate at the parent of J: 

Vc G J,w{J) > ^F(c). 

In addition, Jj -J^fer = Ejep B . t Sj 7$x)> anc ^ ^' therefore, suffices to show that 
for every J G ^Bisection, fj J^fy > 1< Let d G J be such that .F(d) is maximal in J. 
Then 

2cte /" 2di 2 2 F(d) , 



/j " J j F(d) F(d) y ' ~ F(d) 2 

In the case when the Bisection algorithm does not terminate, we can look at the 
partition P at any moment in time. The above argument shows that #P is still 
bounded by the integral Jj2dx/F(x). Since #P can be chosen to be arbitrarily 
large, this shows that the integral is unbounded. □ 

3.2. A stopping function for SqFreeEVAL. The next goal is to transform the 
inequality w(J) < S/ (m(j)) m to a stopping function. Currently, it is not a stopping 
function because the function on the RHS is not for an arbitrary point of J, but 
for a specific point, the midpoint. We begin to turn this into a stopping function 
via the following lemma: 

Lemma 3.1. Let z — (zi, • • • , Zd) with m > and y G M such that y > and 
Zi > y for all i. Then 

(1) HM(z -y)> HM(z) -d-y 

(2) HM(z) < d ■ z, L Mi. 

Proof. For Inequality (flj, we expand each of the harmonic means and get the 
following equivalent inequality: 

d d 

>=^--d-y. 



Noting that all of the denominators are positive, clearing fractions gives that this 
inequality is equivalent to the following inequality: 



y 



yj \ z iJ z i-y 



This inequality is easily justified by combining similar terms on the RHS to obtain 
a sum with general term 

_J 1 = y 

Zi-y z l Zi(zi-y)' 

which is a term that appears on the LHS. Since the remaining terms on the LHS 
are positive, this proves the first inequality. 

For Inequality ([2]) , we expand the harmonic mean to get the following equivalent 
inequality: 

d 
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Once again, the denominator is positive, so by clearing fractions we have that this 
inequality is equivalent to the following inequality: 



Since all of the terms on the RHS are positive and include the term on the LHS, 
this proves the second inequality. □ 

Let Go(x) = 3 y^(x) ' then Go is a stopping function for EVAL: Let J be an 
interval such that J contains x and let m be the midpoint of </; then, \x— m\ < ^ ■ 
Assume now that w{J) < 3^-7^ ■ Then, inequality ([2]) in Lemma |3 . 1 1 implies that 



\x — ot\ > s^z) > ~~^P~ f° r all a € a/. This setup implies the following inequalities: 
(J) < _l MJ) < RM(\x-a f \-^l) < UM(\x-a f \-\x-m\) 1 



2 d d Yif(m) 

The second inequality follows from Lemma 13.11 and the fact that the terms of the 
harmonic mean HM(|a: — a/| — \x — m\) are all positive (because of the bound 
on w(J) above). The remaining inequalities follow from the monotonicity of the 
harmonic mean. The last inequality also uses that \x — ctj\ — \x — m\ < \m — aj\ 
by the triangle inequality. When combined with the observations from Section [2.2[ 
this implies that J is SqFreeEVAL terminal. Similarly, let G\ (x) — 3S where 

£/' is the corresponding function for /', then G\ is also a stopping function for 
the SqFreeEVAL algorithm. Finally, let G(x) — max{Go(x), Gi(x)}, then G is an 
everywhere positive stopping function for the SqFreeEVAL algorithm. 

4. Size of the subdivision tree of the SqFreeEVAL algorithm for the 

benchmark problem 

In this section, we prove that the size of the subdivision tree of the SqFreeEVAL 
algorithm is 0{d(L + \nd)) where L is the number of bits needed to write the 
coe fficients of f. In this case, the absolute value of all the roots is bounded by 
2 L |Ya3, l2000j j (this bound comes from the original /, not from the square free 
substitution). Hence, we can assume wlog that b = —a = 2 L . By Theorem 13. 1[ the 
complexity of the SqFreeEVAL algorithm is bounded by Jj jr^dx. The crossover 
points of G arc difficult to determine, however, so we replace this integral by a 
slightly larger one which is easier to evaluate: For any x £ I, let R x be the set of 
roots in a//' which are closest to x. Similarly, for a € a//', let I a be the set ofxEl 
such that no other root in a//' is closer to x than a. Note that x G la iff ot £ Rx 
and that two of the ia's are either disjoint (except for endpoints) or coincide (in 
the case of complex conjugates). Therefore, these / Q 's determine a partition of /. 
Also, let S be the set of endpoints of the la's; then, for all points x E I \ S, one 
has R x C a.f or R x C ctf because / and /' do not share roots. We define another 
function F(x): 

!G\{x) x S and R x C a/ 
G (x) x<£S and R x C a r . 
G(x) xeS 

Note that although S might not correspond to the crossover points of G, pointwise, 
F(x) < G(x) since G is a maximum of the terms which can occur in F. This implies 
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the following inequalities: 

^ I W) dx - L W) dx - 1 £ —c^, L v~a[ 

For the second inequality let x £ S, then x is either closest to a root of / or a root 
of /'. If x is closest to a root of /, then R x C ay and -p^y = 3S/'(a;). In this case, 
the sum to the right of the inequality includes all of the roots in a f as well as some 
roots in a/. Thus, at least all of the terms of E/'(x) appear on the RHS of the 
inequality. The case where x is closest to a root of /' is similar. This implies the 
inequality because the set of points for which this inequality may fail is a measure 
zero subset of S. 

4.1. Evaluating the integrals. Consider the shape of each of the regions where 
we integrate: since all of the integrals are of the form f* ^fz^y , we evaluate a general 
integral of this form where r and s lie on the same side of Re(a). 
• In the case where a is real: 

if s > r > a if r < s < a 

3dx 



r 3dx r 3dx r m x r 

J r \x- a\ J r x - a J r \x - a\ J r 



a — x 

= 31n(|s - a|) - 31n(|r - a\) = 31n(|r - a|) - 31n(|s - a\ 

These logarithms will be bounded in the next section. 
In the case where a is not real: 



Jr |*^ ^| Jr 

-L 



r ^J{x-Rc{a)) 2 + Im(a) 2 

( S -Re(a))/|Im(a)| 3 



zdx 



(r-Re(a))/|Im(a)| y/y 2 + 1 



dy 



.' s — Re(a)\ . (r — Re(a) 

3 arcsmh — — — , , , — 3 arcsmh ' 



Im(a)| ) \ |Im(a)| 

This is now bounded via the relationship between Re(a) and r, s. If s > 
r > Rc(a), then: 

, ,'s— Re(a)\ „ . , (r — RefaV 
3 arcsmh — — — — — — 3 arcsmh ' 




| Im(a)| y \ | Im(a) 
3 ln(s - Re(a) + y/(s - Rc(a)) 2 + Im(a) 2 ) 



- 3 ln(r - Re(a) + x/(r - Rc(a)) 2 + im(a) 2 ) 
< 31n(2|s-a|) -31n(|r-a|). 
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If r < s < Re(a), then the computation is similar, and the integral is 
bounded above by 31n(2|r — a\) — 31n(|s — a|). These logarithms will also 
be bounded in the next section. 

4.2. Finishing the bound on the SqFreeEVAL algorithm. In this section, we 
use the computation from the previous section to prove the main result of this paper. 
To do this, we consider the roots a £ ctff with two different cases depending on if 
a is real or not. 

• If a is real; then a G I a and let I a = [c, d]. Then, the term corresponding 
to a in the RHS of Inequality © consists of J" ,|% = f 2L j|% + 

Id \ x-a \ ■ Note that integrals may be zero which happens when c = — 2 L 
or d = 2 L . Then, using the bounds derived in the preceding section on 
these integrals, it follows that they are bounded by: 

31n(| - 2 L - a\) - 31n(|c - a\) + 31n(|2 L - a|) - 31n(|d - a\). 

The positive terms are bounded by O(L) (the leading term is 6 ln(2)L) and 
for the negative terms, note that c and d are points which are equidistant 
from a and another root of //', e.g., c is equidistant from a and (3 e oljji 
where Ip is the interval immediately to the left of I a . Then, ln(|c — a|) is 
bounded below by the logarithm of half the distance from a to /3. 

• If a is not real, then the term corresponding a in the RHS of Inequality ([3]) 
consists of Jj.j |-f^y which is bounded above by f T <^_ a i - By splitting this 

integral at Re(a), the integral is equal to J^ { L a) jjj^ + J^ c[a) ||^. Using 
the bounds derived in the preceding section on these integrals, it follows 
that these integrals are bounded by: 

31n(2| - 2 L - a\) - 31n(| Re(a) - a\) + 31n(2|2 L - a\) — 3 ln(| Re(a) - a|). 

The positive terms are bounded by O(L) (the leading term is 61n(2)L) 
and for the negative terms, note that | Re(a) — a\ = | Im(a)|, which is the 
logarithm of half the distance between a and a. 
Combining all of the 0(L)'s which appear in the integrals results in a bound of 
0(dL) (the leading term is 6(ln2(2d— The sum of the logarithmic distances 

between roots are bounded simultaneously vi a the standard Mahler-Davenport 



lower bound on distan ces between roots, see [Davenporti 119851 IDu et al.l . 12007 



Eigenwillig et all 12006] . To do this, we construct a directed graph whose nodes 
are the roots in a ff> and whose edges represent the logarithms which must be cal- 
culated. In this graph, the edges satisfy the conditions of the Mahler-Davenport 
bound and are chosen so that the in-degree of any node is at most 2. For each pair 
of complex roots, (a, a), connect them with two directed edges, one from a to a, 
the other in the opposite direction. On the other hand, if a is real, then let /3 be 
a root where Ip lies immediately to the right of I a and 7 be a root where 1 1 lies 
immediately to the left of I a (provided I a is not the rightmost or leftmost interval 
in the partition, respectively). Those of /3 or 7 which are real are connected to a 
so that the arrow points in the direction of decreasing absolute value. If /3 is not 
real, then connect a to either ft or /3, whichever has positive imaginary component. 
On the other hand, if 7 is not real, then connect a to either 7 or 7, whichever has 
negative imaginary component. Again, these edges are directed so that the arrow 
points in the direction of decreasing absolute value. By inspection, we find that the 
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maximum in-degree of this directed graph is 2. The Mahler-Davenport bound can 
then be applied twice to find the result. The bound implies that the sum of the 
negative logarithmic distances between the roots appearing in this construction is 
bounded above by: 

12 • In ( 1 M(ff) 2d - 2 ( ^ 1 (2d - l) d ) . 

The discriminant will be an integer and therefore the discriminant term is bounded 
above by 1. The Mahler measure of //' is bounded in terms of the 2- norms of 
the coefficients of the original / and /': M(ff) = M(f)M(f) < ||/|| 2 ||/'||2 < 
(2 L y/d + l)(d2 L \fd). Therefore, this portion is bounded by 0(dL + d\nd) (the lead- 
ing term is bounded 241n(2)dL+42dlnd)). Thus, the complexity of the SqFreeEVAL 
algorithm is 0(d(L + hid)) (the leading term is bounded 361n(2)cLL + 42c£lnd < 
25dL + A2d\nd). 

If / or /' was replaced by a square free version, we used the original / and 
/' because the square free versions of / and /' divide the original functions, and, 
therefore, the Mahler measure of the product of the square free versions is bounded 
above by M(ff'). In fact, the 2-norms of the coefficients of the original functions 
are often smaller than the 2-norms of the square free versions. 



5. Conclusion 



In this paper, we provided a complexity analysis of the SqFreeEVAL algorithm 
and showed it to be optimal under the weak assumption that L > In d. To ac- 
complish this, we used the novel technique of continuous amortization through 
stopping functions. The simplicity of this argument exhibits the utility of this 
techniqu e: the proof of the next c losest complexity bound for an EVAL-type algo- 
rithm in SagralofF and Yarll2009j is significantly more co mplex 



1966, Mitchell 



The SqFreeEVAL algorithm is very easy to implement Moord . 
199(1 iPlantinga and Ve gtcr , I2Q04L IPlantingal . 120061 iKamathl l2010j and it now joins 
the Sturm and Descartes methods by having a subdivision tree which grows at the 
rate 0(d(L+lnd)). It, therefore, may become more prevalent in practical situations 
because it has several desirable properties. This also answers a question raised in 
Henricil [l970j concerning the good behavior of this technique. 



In addition, the continuous amortization technique can be used to bound the 
number of subdivisions over any interval, and, therefore, may find many more 
applications for different types of questions about subdivision algorithms. For ex- 
ample, in many practical applications, the question is to find the roots in a given 
domain, not just for the benchmark domain; continuous amortization may provide 
a comparison of different algorithms in these situations. 

We close with some continuing research and questions: 

• The algorithm for finding complex roots appearing in Sagraloff and Yapl 
2009] is very similar to the SqFreeEVAL algorithm. We are currently prepar- 



ing a simplification of their work using the results from this paper. 
There are many bisectio n algorithms where continuous amortization may be 
useful, see, for example, Henricil. 1970 L Yakoubsohn L 2005, Sagral off and Yai 



Plantinga and VegterL I2004L IPlantingaL 120061 ISnvderl ll992L iGalehous- 



2009 



2009 



Burr et al.Ll2010llEigenwillig et a l.. 2006, Du et all . 120071 



We plan on extending our techniques to these cases. In particular, 
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stopping f unctions which are appropriate for the two dimensional cases 
treated in [Plantinga and Vegterl . I200I IPlantingaL l2006l iGalehousd . l2009j | 
would be very useful because current techniques have not been fruitful in 
establishing complexity bounds of these algorithms. 

If /' was not square free, then the test for condition (C\) in Algorithm 2.1 
is based on the square free part of /', not the original function. The 
SqFreeEVAL algorithm, however, will continue to terminate and be correct 
even when this substitution does not occur, i.e., when the original /' is used. 
For this reason, it is likely that the above substitution is extraneous. For 
example, in the simplest cases where /' is not square free and the integral 
in Inequality ([3]) can be calculated by hand, the result is 0(d(L + hid)). 
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